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Abstract 

We sequenced mitochondrial genes ( COI, COII, Cyt-b ) of accepted Latin America 
tapir species {Tapirus pinchaque, T. terrestris and T. bairdii) as well as an alleged new 
species, T. kabomani. 

The mountain tapir (T. pinchaque) is a relatively rare large mammal species. Some 
population censuses indicate that no more than 2,000 mountain tapirs are left in the 
wilderness areas of Colombia, Ecuador and Peru. Our results showed that the gene 
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diversity levels are medium to low with respect to other mammals sequenced for the 
same or similar genes. However, these gene diversity levels are not impoverished, which 
means that the genetic situation of this species is not as critical as its population censuses 
suggest. It will be crucial to determine the gene diversity levels in certain populations not 
included in the current work (the eastern and, possibly, western Andean Cordilleras in 
Colombia as well as the Tabaconas Namballe National Sanctuary in Peru), because they 
are probably the smallest populations of this species. On the other hand, the lowland tapir 
(T. terrestris ), the species with the largest geographical distribution in Latin America, 
showed the highest gene diversity levels of all the other tapir species studied. 
Additionally, the genetic structure of T. terrestris is clearly more robust than that of T. 
pinchaque. Different geographic populations of both species showed different 
demographic trends throughout time. Our results including five samples of T. kabomani 
showed this taxon to be a haplogroup within T. terrestris , reducing the likelihood of T. 
kabomani being a new full species. Finally, we also analyzed the influence of diverse 
Pleistocene climatic changes on the mitochondrial haplotype diversification of T. 
terrestris and T. pinchaque. The Pleistocene Refugia and the Recent Lake hypotheses 
probably played integral roles in the evolutionary history of T. terrestris. In contrast, the 
Pleistocene Refugia hypothesis involving the Andes, which probably played an important 
part in the genetic diversification of other mammals, did not have a significant impact on 
T. pinchaque. 

Keywords: Mitochondrial genes (CO!, COII and Cyt-b), Tapirus, T. pinchaque, T. terrestris, 
T. bairdii, “T. kabomani ” conservation genetics, speciation. Pleistocene biodiversity 
hypotheses 


Introduction 

Large mammals are often considered umbrella or emblematic species whose presence 
may greatly affect community structure (Meffe and Carrol, 1997). In general, the extraction 
of large vertebrates, especially those at the top of the trophic chain, can provoke strong 
extinction waves in complex natural systems (Pimm, 1991). These species can be considered 
as key to overall community structure because they live sympatrically with many other 
species. For example, East (1981) showed that the preservation of viable populations of the 
African wild dog ( Lycaon pictus) and the cheetah (Acinonyx jubatus), due to their large 
hunting territories, helped to preserve the communities of other large African mammals. 

In Latin America, three wild large perissodactyla species have been traditionally 
accepted: the Baird or Central America tapir ( Tapirus bairdii), the mountain tapir (T. 
pinchaque ) and the lowland or Brazilian tapir ( T. terrestris). 

The fossil evolution of the Tapiridae is old and very interesting. Perissodactyla is a very 
old group of mammals that arose around 60 Million years ago (MYA). In the fossil record, 
there are representative specimens from five main superfamilies (Tapiroidea, 
Rhinocerotoidea, Chalicotheroidea, Equoidea and Brontotheroidea) and 14 different families 
(Savage and Long, 1986; Holbrook, 1999). The order reached its maximum diversity peak 
during the Eocene and then declined from 14 to 4 families during the upper Oligocene 
(Radinsky, 1969; Froehlich, 1999; MacFadden, 1992; Metais et al., 2006). The Tapiridae 
family (Gray 1821) is currently composed of one unique genus, Tapirus (Brisson 1762). 
Colbert (2005) defined the Tapiridae family as the clade conformed by the most recent 
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common ancestor of Protapirus. The oldest record of Tapirus comes from the European 
Oligocene (33-37 MYA), where the fossil remains are found until the Pleistocene (McKenna 
and Bell, 1997). In North America, the Tapirus records indicate that they were present in the 
Middle Miocene through the present (Hulbert, 1995), whereas for Asia the records indicate 
that Tapirus has been in existence since the Lower Miocene (Deng, 2006). Currently, 
additional to the three quoted Latin American tapir species, there is a fourth species in 
Southern Eastern Asia, the Malayan tapir (T. indicus). Two of these species are exclusively 
distributed in South America and they are the subject of this work, from a genetics point of 
view. The mountain tapir (T. pinchaque) has a unique and critical role because it is an 
indispensable seed disperser in the Andean mountains (Colombia, Ecuador, Northern Peru) 
but its conservation status is in jeopardy because of its relatively low numbers (less than 
2,500 individuals; Cavelier et al., 2010). Additionally, this is an umbrella species with an 
average territory size per tapir of 880 ha (Downer, 1996). This is a substantially large home 
range. Downer’s work in the Sangay National Park suggested that a minimum of 1,000 
mountain tapirs were needed to ensure a reproductively viable population. Such a population 
would need at least 293,500 ha. The necessary protection of this area would positively benefit 
many other sympatric animal species, like the Andean bear, puma, and various deer species 
(Downer, 1996). Indeed, the cloud forests of the Andean region, where the mountain tapir 
lives, contain some of the greatest organismic diversity on Earth (Brehm et al., 2005; Kromer 
et al., 2005), with a noteworthy number of endemic species (Kessler, 2002). However, about 
90 % of the mountain Andean forests have been deforested (Henderson et al., 1991) and a 
great proportion of paramos has been transformed via burns into extensive ranchlands and 
croplands (potatoes) (Verweij, 1995). Additionally, around 83 % of Colombia’s mountain 
forests have been highly affected by human activity (Cavelier and Etter, 1995). Besides 
ranching and agriculture, the mountain tapir is also highly threatened by poaching and illegal 
trade of its body parts (Downer, 2003). In Andean areas, the introduction of livestock has also 
spread new diseases and attracted more potential predators of the mountain tapir. 

The lowland tapir has a wide geographical distribution across a great part of South 
America. From an ecological perspective, the lowland tapir is considered the “architect” of 
their Neotropical habitats because it is a very efficient disperser and seed predator (Wallace et 
al., 2010). Cabrera (1961) considered four morphological and geographical subspecies of T. 
terrestris : T. t. colombianus, T. t. terrestris, T. t. aenigmaticus, and T. t. spegazzinni. T. t. 
colombianus (Hershkovitz 1954; Type locality: El Salado, between Valencia and Pueblo 
Bello on the eastern slope of the Sierra Nevada de Santa Marta, in Magdalena, Colombia) 
inhabits Northern Colombia and around Maracaibo Lake in Venezuela. T. t. terrestris 
(Linnaeus 1758; Type locality: Pernanbuco, Brazil) lives in Venezuela, Guyana, Suriname, 
French Guyana, with its major fraction in Brazil (including Amazonas) up to the Misiones 
within Argentina. T. t. aenigmaticus (Gray 1872; Type locality: Macas, Eastern Ecuador) is 
distributed in Southeastern Colombia, Eastern Ecuador and Northern/Eastern Peru 
(Amazonian area). T. t. spegazzinni (Ameghino 1909; Type locality: Rio Pescado, Oran, 
Salta, Argentina) is distributed in Southeastern Brazil and Mato Grosso, Paraguay, Eastern 
Bolivia and Northwestern Argentina. 

Recently, Cozzuol et al., (2013) claimed the existence of a new tapir species (“T. 
kabomanr) in South America based on a comparison of mitochondrial Cytochrome-b ( Cyt-b) 
gene sequences of four Amazon tapir individuals (two Brazilian animals sampled by these 
authors and two Colombian specimens sampled by the first author of the present study, M. R- 
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G) to the 45 Cyt-b gene sequences published by Thoisy et al., (2010). They also enclosed the 
results of two additional mitochondrial genes (mitochondrial Cytochrome Oxidase subunit I, 
COI, and subunit II, COIT) for six T. terrestris, one T. pinchaque and three individuals of the 
alleged new species “T. kabomani." They concluded that the three mitochondrial genes 
supported “T. kabomani” as a full species. Voss et al., (2014), however, questioned the 
validity of “T. kabomani" as a different species from T. terrestris. 

The current work provides some insight about the genetic structure and heterogeneity, 
demographic evolution, spatial structure, biological conservation, and systematics of the Latin 
American tapirs. Moreover, this work also explores the influence of the Pleistocene period on 
the Latin American tapirs. And, all of these insights are based on the analysis and study of 
mitochondrial genes Cyt-b, COI, and COII. These new data should be helpful in validating or 
rejecting the claim of a new species, “T. kabomani". 


Material and Methods 

Two sets of data were used in this study. The first set consisted of 93 Latin American 
tapir individuals analyzed for three specific concatenated mitochondrial genes ( Cyt-b + COI + 
COII). They were the same genes analyzed by Cozzuol et al. (2013) and whose mitochondrial 
sequences are accessible in GenBank for the two Brazilian specimens they classified as an 
alleged new species “T. kabomani .” Breaking the first group down by taxa, it had 46 T. 
pinchaque, 29 T. terrestris, 13 T. bairdii, and 5 T. kabomani. T. pinchaque was constituted by 
26 from Ecuador representing three populations, Northern Ecuador, Sangay National Park 
and Podocarpus National Park, and the remainder of the 46 (20) was from Colombia 
representing three populations: Los Nevados National Park, Tolima and Purace National 
Park. 

Another 29 of the first set were T. terrestris. Five of these individuals were from 
Ecuador, four were from Colombia, seven from Peru, three from Bolivia, three from Brazil, 
three from Surinam, and one was from Argentina. Three animals of the 29 were from 
American zoos in Cincinnatti, Milwakee and Columbus. Another 13 of the 93 were T. bairdii. 
Nine of these animals were from Panama, three were from Costa Rica and one was from 
Colombia. Five of the 93 were T. kabomani- including the two Brazilian specimens of 
Cozzuol et al., (2013).We detected three additional specimens with mitochondrial haplotypes 
of the alleged T. kabomani. There were from 1) San Martin de Amacayacu along the Amazon 
River in Colombia, 2) the Mazan River which is a tributary of the Napo River, in the Peruvian 
Amazon, and 3) near Tena along the upper Napo River in the Ecuadorian Amazon. 

However, at least, two of these animals (the Colombian and the Ecuadorian ones) 
presented typical morphotypes of T. terrestris, although they showed “T. kabomani ” 
haplotypes (see Figure 1). The Peruvian sample consisted of skin obtained from hunters and 
we could not determine the morphotype of the specimen. 
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Figure 1. An individual with morphotype of Tapirus terrestris but with “T. kaboman” mitochondrial 
haplotype from San Martin de Amacayacu, Amazon River, Colombian Amazon (A); An individual 
with morphotype of Tapirus terrestris but with “T. kaboman ” mitochondrial haplotype from near to 
Tena, upper Napo River, Ecuadorian Amazon (B). 

The second set of individuals contained 141 T. terrestris, sequenced at the Cyt-b gene. 
Forty-one of these came from different regions of Colombia. Of these, 1 was from Bajo Sinu- 
Tierra Alta, (Cordoba Department), 2 were from Mesay River (Caqueta Department), and 
another was from Fondo Canaima, (Vichada Department). Of the 41, 18 were also from 
Leticia to San Juan de Atacuarf (Amazonas Department), 7 from the Eastern Colombian 
Llanos (Meta Department), 3 from Pto. Inirida (Guania Department), 3 from Palomino River- 
Sierra Nevada de Santa Marta and 6 from the Antioquia Department. Five of the 141 were 
from Venezuela (El Zulia, Maracaibo). Eleven of the 141 were from French Guiana (Camopi 
River). Seven animals were from Ecuador (4 from Limoncocha, Sucumbios and 3 from Coca, 
Sucumbios). Thirty animals were from Peru [(one from Arica (Curaray River), 2 from Napo 
River (Nueva Vida and Mazan), 7 from Nanay River, one from Requena (Ucayali River), one 
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from Bretana (Canal del Puhinauva-Ucayali River), 4 from Pucallpa (Ucayali River), and 15 
from Pto. Maldonado (Madre de Dios River)]. Eleven animals were from Bolivia (9 from 
Mamore River, one from Chimore River and one from Villa Bella at the Beni River). Twenty 
four animals were from Brazil [2 from Yavarf River, 12 from Tabatinga, 2 from Negro River, 
2 from Santarem (Para state) and 6 from the Amazon’s mouth (Para state)]. One animal was 
from Paraguay (from Hemandarias). And, 4 animals were from Argentina (the Yungas area in 
Salta-Jujuy). Additionally, one animal was from the Barcelona Zoo (Spain) and 5 animals 
were from US zoos. One animal of unknown origin was also analyzed. 


Molecular Procedures 

DNA from teeth, bones, muscle, and skins, were obtained with the phenol-chloroform 
procedure (Sambrook et al., 1989), whereas DNA samples from hair and blood were obtained 
with 10% Chelex® 100 resin (Walsh et al., 1991). Amplifications for Cyt-b gene were 
achieved using primers L7 (5' ACC AAT GAC ATG AAA AAT CAT CGT T 3') and H6 (5' 
TCT CCA TTT CTG GTT TAC AAG AC 3'), that had been designed for perissodactyles 
(Tougard et al., 2001). The PCR reactions were performed in a 50-pi volume, including 10 pi 
of lOx Buffer, 7 pi of 3 mM MgCl 2 , 2 pi of 10 mM dNTPs (dNTP Mix Promega), 4 pi (15 
pmol) of each primer, one unit of Taq DNA polymerase (genTaq polimerasa), 2 pi of DNA 
from blood, skin or muscle tissue (50-200 ng/pl) or 4-10 pi of DNA from hair and teeth (6-35 
ng/pl) and a variable quantity of ddH 2 0. PCR reactions were carried out in a Geneamp PCR 
system 9600 (Perkin Elmer) and in an iCyclerTM BioRad thermocycler. We used the 
following temperatures: 94 °C for 5 minutes, 35 cycles of 50 s at 94 °C, 50 s at 53 °C and 1.5 
minutes at 72 °C and a final extension time for 10 minutes at 72 °C. For the COI and COll 
genes, the amplification conditions followed Ashley et al., (1996) and Hebert (2003, 2004). 
All amplifications, including positive and negative controls, were checked in 2 % agarose 
gels, employing the molecular weight marker 4>X174 DNA digested with Hind III, Hinf I and 
HyperLadder IV. The gels were visualized in a Hoefer UV Transilluminator. Those samples 
that amplified were purified using membrane-binding spin columns (Qiagen). The double- 
stranded DNA was directly sequenced in a 377A (ABI) automated DNA sequencer. The 
samples were sequenced in both directions and all the samples were repeated to ensure 
sequence accuracy. 


Data Analyses 

Genetic Diversity, Linkage and Heterogeneity Analyses 

The sequences were edited and aligned with BioEdit Sequence Alignment Editor (Hall, 
2004) and DNA Alignment (Fluxus Technology Ltd). 

We used the following statistics to determine the genetic diversity at the three 
concatenated mitochondrial genes for T. terrestris, T. pinchaque, T. bairdii and the alleged 
“T. kabomani the number of polymorphic sites (S), the number of haplotypes (H), the 
haplotypic diversity (H d ), the nucleotide diversity (n), the average number of nucleotide 
differences (k) and the 9 statistic by sequence. 
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The possible linkage disequilibrium within the Cyt-b gene for T. terrestris was evaluated 
by two different statistics. We obtained the average value of R~ for all the comparison pairs 
among nucleotide sites (ZnS) (Kelly, 1997) as well as among the adjacent polymorphic 
nucleotide sites (Za), and ZZ (Rozas et al., 2001). This is equal to Za - ZnS, using the 
software DNAsp 5.10 (Librado and Rozas, 2009). 

Different tests were carried out to measure genetic heterogeneity, and possible gene flow 
estimates, among the four Latin American tapir taxa studied. These tests were those of 
Hudson et al., (1992a,b) (H S t , K S t , K S t*, Z, Z*), Hudson (2000)’s Snn test and the chi- 
square test on the haplotypic frequencies with permutation tests using 10,000 replicates. We 
also estimated the Gst statistic from the haplotypic frequencies and the yst, Nst and F ST 
statistics (Hudson et al., 1992a) from the nucleotide sequences. 

We carried out two different AMOVA analyses at the Cyt-b gene for T. terrestris to 
determine if the gene diversity within this species was distributed in different hierarchical 
geographical levels (Excoffier et al., 1992). The first analysis was carried out taking into 
account the six haplogroups found in the phylogenetic analysis clustered in three main 
population sets: northern, Amazonian, and southern. Thus, this analysis was applied to the 
overall geographical T. terrestris range analyzed. The second analysis was completed taking 
into account six different geographical areas within the Amazon basin: 1-Northwestern 
Amazon in Colombia and Brazil, 2- Central and Eastern Brazilian Amazon, 3-Napo River and 
a tributary in Peru and Ecuador, 4- Ucayali River in Peru, 5-Madre de Dios River in Peru and 
6- Mamore River and a tributary in Bolivia. These six geographical areas were clustered into 
two main groups: northern Amazon (the first three areas) and southern Amazon (the last three 
areas). This analysis was carried out to determine the gene diversity structure in the Amazon 
region. The fixation indices of Wright (1951) were estimated: ® sc (variation of populations 
within the groups), <D c t (variation among groups) and ® st (variation among individuals). 
These analyses were carried out by means of the software ARLEQUIN 3.0 (Excoffier et al., 
2005). 

Phylogenetic Analyses 

The Modeltest (Posada and Crandall, 1998) and the Mega 5.1 software (Tamura et al., 
2011) were applied to determine the best evolutionary nucleotide model for the analyzed 
concatenated gene sequences for all the Tapirus taxa studied. 

To determine if T. pinchaque or T. terrestris, contained a more conspicuous genetic 
structure, we used a consensus maximum parsimony (MP) tree for the three concatenated 
mitochondrial genes. We also obtained a maximum likelihood (ML) tree to determine how 
many significant clades were within T. terrestris, at the Cyt-b gene. In this analysis, 80 T. 
terrestris haplotypes were considered (including two haplotypes later classified as “T. 
kabomani”). To carry out this task, we tested the hypothesis that the T. terrestris sequences 
fall into one to 10 different groups. These classifications were contrasted with the ML tree. 
For this, we performed parametric bootstrapping and a posteriori significance test with the 
Swofford-Olsen-Waddell-Hillis test (SOWH; Huelsenbeck and Bull, 1996; Swofford et al., 
1996). The 10 hypotheses were used as a model tree for parameter estimation and for 
generating 100 replicate data sets in the software Seq-Gen 1.2.5 (Rambaut and Grassly, 1997) 
which presented a uniform base composition. Goldman et al., (2000) demonstrated that this 
procedure can increase power in rejecting the null hypothesis and is better than typical 
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nonparametric tests for comparisons of a posteriori hypotheses. The same was performed with 
the Shimodaira and Hasegawa (1999) test (a nonparametric SH test). 

We conducted two additional phylogenetic trees. One was a ML tree with the 93 tapirs 
sequenced at the three concatenated mitochondrial genes to determine the relationship of the 
five “71 kabomani” with regard to the other tapir taxa. The second was an ML tree with some 
animals only sequenced at the COI gene since this gene is used as a barcode to discriminate 
different species (Herbert et at., 2003, 2004). Similarly, we estimated the Kimura 2P genetic 
distance (Kimura, 1980), for the three mitochondrial genes, among all the Neotropical 
Tapirus taxa because this genetic distance is a standard measurement for barcoding tasks 
(Hebert et al., 2003, 2004). 

Possible divergence times were obtained using a Median Joining Network (MJN) 
(Bandelt et al., 1999). These were applied to all the haplotypes of the three concatenated 
mitochondrial genes for all the tapir taxa considered as well as to the T. terrestris haplotypes 
at the Cyt-b gene. These were constructed with Network 4.6 software (Fluxus Technology 
Ltd). The p statistic (Morral et al., 1994) was estimated and transformed into years. The 
standard deviation of p was also calculated (Saillard et al., 2000), which is unbiased and 
highly independent of past demographic events. Thoisy et al., (2010) used two different 
mutation rates: 5.6 x 10 and 2.5 x 10 substitutions/site/million years respectively. Relative 
to our work with the genus Tapirus, the first rate translated into around one mutation each 
537,634 years. The second rate was equal to about one mutation each 120,482 years. This last 
mutation rate agrees quite well with that reported by Nabholz et al., (2008, 2009) for 
mammals at the Cyt-b gene. We employed an intermediate value of 1.5 x 10~" 
substitutions/site/million years (around one mutation each 200,000 years), closer to the 
second value used by Thoisy et al., (2010). 

Demographic Changes 

We used two methods to determine possible demographic changes across the natural 
history of T. pinchaque, T. terrestris and “T. kabomani ” for the three concatenated 
mitochondrial genes as well as for the Cyt-b gene in the lowland tapir. (1) Following the 
method of Rogers and Harpending (1992) and Rogers et al., (1996) we used a mismatch 
distribution (pairwise sequence differences). The raggedness rg statistic (Harpending et al., 
1993; Harpending, 1994) was used to determine the similarity between the observed and the 
theoretical curves. (2) We used the Fu and Li D and F tests (Fu and Li, 1993), the Fu F s 
statistic (Fu, 1997), the Tajima D test (Tajima, 1989) and the R 2 statistic (Ramos-Onsins and 
Rozas, 2002) to determine possible changes in population size (Simonsen et al., 1995; 
Ramos-Onsins and Rozas, 2002). All of these statistics and tests were obtained by means of 
the DNAsp 5.10 and Arlequin 3.0 programs. 

Spatial Structure 

Garnering this information can assists in our understanding of the evolutionary events 
that have determined the natural history of the tapir species. A Mantel’s test (Mantel, 1967) 
was used to detect possible overall relationships between a genetic matrix among individuals 
(Log-Det genetic distance; Nei and Kumar, 2000) and the geographic distance matrix among 
the individuals analyzed. In this study, Mantel’s statistic was normalized according to 
Smouse et al., (1986). This procedure transforms the statistic into a correlation coefficient. 
The geographic distances were measured with the Spuhler’s (1972) procedure, where D = 
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arcos (cos X (i) . cos X ( j, + sin X (i) . sin X^) cos IY (i ) - Y (JI I), where X (n) and Y (n) are the latitude 
and longitude of the nth individual sampled, respectively. The significance of the correlations 
obtained was tested using a Monte Carlo simulation with 1,000 permutations. This test was 
undertaken by means of the software NTSYS v. 2.1 (Rohlf, 2000). This analysis was carried 
out for T. pinchaque, T. bairdii and the overall sample of T. terrestris for the three 
concatenated mitochondrial genes as well as for the different haplogroups found and for 
different geographical regions at the Cyt-b gene for T. terrestris. For all the Amazon T. 
terrestris individuals sequenced at the Cyt-b gene, an AIDA analysis was used (Bertorelle and 
Barbujani, 1995). The expressions of the respective AID As coefficients are as follows: 

— (it ^ i =1 to n ^ j > i to n Wij E k = j t 0 s (Pik Pk) (Pjk Pk)) / 

(W Y i =1 to n ^ k =1 to S (Pik — Pk) ) 


and 

2 

CC — ((ft " l)^i=lton^j>iton Wjj Y k =1 to S (Pik — Pjk) ) / 

(2W E i = 1 to n £ k =1 to S (Pik - Pk)”) 

where n is the sample size, W is the number of pairwise comparisons of a distance class 
given, p ik and p jk are the haplotypes of the ith and jth individuals, respectively. At the kth 
nucleotide site, p k is the kth element of the average vector and w,, is one if individuals i and j 
are within the same distance class; otherwise it is zero. Summation includes the E nucleotide 
sites for all the n individuals analyzed. To connect the individuals sequenced within each 
distance class, the Gabriel-Sokal network (Gabriel and Sokal, 1969; Matula and Sokal, 1980) 
and the Delaunay's triangulation with elimination of the crossing edges (Ripley, 1981; Upton 
and Fingleton, 1985; Isaaks and Srivastava, 1989) were used. However, the results were very 
similar in each case. The Bonferroni (Oden, 1984), Oden's Q and the Kooijman’s tests were 
estimated to determine the statistical significance of the autocorrelation coefficients. 


Results 

Gene Diversity, Genetic Heterogeneity and Phylogenetic Considerations for 
the Latin American Tapir Taxa 

The GTR with gamma distributed rate variation among sites for both the maximum 
likelihood and the AIC criteria was the best fit nucleotide substitution model at the three 
concatenated genes. 

Out of 4,753 comparisons among the total polymorphic sites at the Cyt-b gene in T. 
terrestris, only 303 comparisons were statistically significant by means of a Fisher exact test 
(6.4 %), but they were not different from the Type I error of 5 %. Similarly, the values of the 
statistics of Kelly (1997) (ZnS = 0.0219; 95 % of Confidence interval. Cl = 0.01498 - 
0.31650), of Rozas et al., (2001) (Za = 0.0244; Cl = 0.01706-0.32607) and ZZ = 0.0026 (Cl 
= -0.05363 - 0.05625) also showed that there was no significant association among the 
polymorphic sites at this mitochondrial gene. Therefore, no evidence of linkage among 
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polymorphic sites was detected at this mitochondrial gene. Out the four taxa of Tapirus 
considered at the three concatenated mitochondrial genes, clearly 71 terrestris showed the 
highest levels of gene diversity (H d = 0.985; n = 0.0094; k = 22.98), followed by 71 
pinchaque (H d = 0.957; n = 0.0041; k = 9.98). In contrast, “71 kabomani” (Hd = 0.900; n = 
0.0033; k = 8.00) and 71 bairdii (H d = 0.910; n = 0.0017; k = 4.08) yielded the lowest levels 
of gene diversity (Table 1). The “71 kabomani ” gene diversity should be enclosed within the 
gene diversity of 71 terrestris (proof 1). The genetic heterogeneity among these four taxa was 
highly significant and of a large magnitude for all the statistics used (for instance, y ST = 0.808 
and F st = 0.903; Table 2a), with the associated indirect gene flow values practically non¬ 
existent (Nm = 0.05-0.12). 


Table 1. Gene diversity statistics for Tapirus terrestris, T. pinchaque, T. bairdii and the 
alleged new species “71 kabomani ” at the mitochondrial genes sequenced ( Cyt-b + COI + 
COII). The statistics estimated were the number of haplotypes (NH), the haplotypic 
diversity (H d ), the nucleotide diversity (n), the average number of nucleotide differences 
(K) and the 0 statistic (= 2N e p; N e = effective female population size; p = mutation rate 

per generation) by sequence 



NH 

H d 

n 

K 

9 per 
sequence 

Tapirus terrestris 

24 

0.985 

±0.014 

0.0094 

±0.0007 

22.985 

±10.42 

29.538 

±9.448 

Tapirus pinchaque 

31 

0.957 

±0.021 

0.0041 

±0.0005 

9.981 

±4.648 

13.652 

±4.161 

“Tapirus kabomani” 

4 

0.900 

±0.161 

0.0033 

±0.0009 

8.000 

±4.483 

7.680 

±4.165 

Tapirus bairdii 

9 

0.910 

±0.068 

0.0017 

±0.0003 

4.077 

±2.173 

5.156 

±2.268 


Table 2. Genetic heterogeneity statistics applied to Latin American tapirs: among four 
Latin American tapir taxa (71 pinchaque, T. terrestris, T. bairdii and the alleged new 
species “71 kabomanV’) (A); between T. pinchaque and T. terrestris (B). *Significant 

Probability (P < 0.01) 

A 

Among four Latin American tapir taxa 


Genetic differentiation estimated 

p 

Gene flow | 

r= 279.000 df = 201 

0.0002* 



H st = 0.0287 

0.0000* 

Yst — 0.8079 

Nm = 0.12 

K st = 0.8026 

0.0000* 

N ST = 0.9078 

Nm = 0.05 

K st * = 0.3433 

0.0000* 

F st = 0.9034 

Nm = 0.05 

Z s = 941.8106 

0.0000* 



Z s * = 6.4554 

0.0000* 


S m = 1.0000 

0.0000* 



B 
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T. pinchaque sample vs. T. terrestris sample 


Genetic differentiation estimated 

P 

Gene flow | 

X 2 = 104.041 df = 70 

0.0052* 

G sx = 0.0411 

Nm= 11.68 

H ST = 0.0414 

0.0004* 

Yst = 0.4389 

Nm = 0.64 

K ST = 0.4167 

0.0000* 

N ST = 0.5158 

Nm = 0.47 

K st * = 0.2421 

0.0000* 

F st = 0.5145 

Nm = 0.47 

Z s = 678.8785 

0.0000* 



Z s * = 6.1920 

0.0000* 


S nn = 0.7876 

0.0000* 



In another heterogeneity analysis, T. pinchaque was only analyzed with regard to T. 
terrestris (Table 2b). In this case, the seven genetic heterogeneity tests were also highly 
significant at the P < 0.05 level. The relative genetic differentiation tests yielded elevated 
amounts of genetic heterogeneity between the Tapirus taxa (yst = 0.439 and Fst = 0.514). The 
gene flow estimates were clearly lower than 1 (Nm = 0.47-0.64), which showed that these 
taxa are genetically disconnected. Although, “T. kabomani ” showed significant differences 
with T. terrestris as well as with T. pinchaque, the genetic differences with T. terrestris are 
considerably lower than with T. pinchaque (six significant heterogeneity tests, ysT = 0.151 
and Fst = 0.454, Nm = 2.82-0.60 versus seven significant heterogeneity tests, ysT = 0.300 and 
F S t = 0.711, Nm = 1.16-0.20, respectively). Therefore, “T. kabomani ” is less differentiated 
genetically speaking from T. terrestris than from T. pinchaque and it showed less 
differentiation in regard to T. terrestris than the genetic differentiation between T. terrestris 
and T. pinchaque. This provides strong evidence (proof 2) that “T. kabomani ’ is more closely 
related to T. terrestris than to the other tapir species, which disagrees with the point of view 
of Cozzuol et al., (2013). 

The genetic structure within T. terrestris is much more robust (greater) than within T. 
pinchaque. This is verified by the consensus MP tree (Figure 2). T. pinchaque showed some 
small geographical clusters with several individuals of nearby geographical regions. Within T. 
pinchaque there were five Colombian individuals (four from Los Nevados National Park and 
one from Gaitania, Tolima) and two individuals from the Alto Papallacta River in the Napo 
Province in Ecuador. It also contained seven Ecuadorian individuals (six from the Napo 
Province and one from Sangay National Park), two from the Napo Province and an 
Ecuadorian cluster of six. 

The Ecuadorian cluster represented three different Ecuadorian populations. Two T. 
pinchaque sequences were more related to those from T. terrestris than to those of the other 
mountain tapirs. These were the cases of one Colombian individual from La Planada, Tolima 
(Colombia) and one Ecuadorian individual from Chaco, Oyacachi, Napo. However, the 
relationships among the T. pinchaque were considerably less structured than those observed 
within T. terrestris. This second species was composed by more robust haplogroups than 
those detected in the mountain tapir. Even an exact probability test with a Markov chain 
length of 100,000 steps only revealed one significant population comparison pair: Los 
Nevados vs. Sangay National Park (p = 0.0297 + 0.0013) out of the six populations of T. 
pinchaque studied. This may indicate that (proof 3) T. terrestris is an older species than T. 
pinchaque. 
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T. terrestris 
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Figure 2. Consensus maximum parsimony tree with 45 Tapirus pinchaque, plus 17 T. terrestris by 
using the concatenated sequences of three mitochondrial genes ( Cyt-b + COI + COIT). 

The most developed genetic structure in T. terrestris can be seen in the AMOVA 
analyses carried out. The first AMOVA, with the six haplogroups, showed that the major part 
of the genetic variance was among lineages (64%; ® sc = 0.762, P-value = 0.000+0.00), 
followed by the gene variance among the individuals (36%; ® st = 0.641, P-value = 0.000 + 
0.000). In contrast, the genetic variance among the three main groups (North, Amazonian and 
South) did not show significant variance (0%, ® ct = -0.510, P-value = 1.000 + 0.000). 
However, in the second AMOVA, by geographical areas within the Amazon basin, the major 
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part of the genetic variance was explained by the difference among the individuals (90%, ® st 
= 0.499, P-value = 0.000 + 0.000). The genetic variance among the populations, although 
significant, only explained 10 % of the genetic variance (® sc = 0.148, P-value = 0.000 + 
0.000). The gene variance between the two main groups was again not significant (0%, O ct = - 
0.057, P-value = 0.913 ± 0.008). 
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Figure 3. Maximum likelihood (ML) tree with the 80 haplotypes (from haplotype 2 to haplotype 81) 
found at the Cyt-b gene (1,140 bp) for the 141 lowland tapirs sequenced analyzed. Numbers in the 
nodes are bootstrap percentages. 
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The lower gene variance among geographical regions relative to among haplogroups is 
explained because the diverse haplogroups coexisted sympatrically at the same points of the 
Northern Peruvian Amazon (four out of six different haplogroups). This was also true for the 
Colombian Amazon (six of six). Northwestern Brazilian Amazon (three out of six), Southern 
Peruvian Amazon (three out of six) and Bolivian Amazon (three out of six). The existence of 
extremely, well-differentiated haplogroups within 71 terrestris could help us to understand 
why “71 kabomani” is an additional haplogroup within “71 kabomani” rather than a full 
species (proof 4). An ML tree, applied to 71 terrestris at the Cyt-b, showed that six 
haplogroups were detected within this species (Figure 3). The SOWH tests (and also the SH 
tests) indicated that there was neither support for the taxonomic schemes of 1 to 5 clusters nor 
for the 7 to 10 clusters. The maximum likelihood trees were significantly different at the 
0.001 level (since 12,569 to 689,999 log likelihood units). However, the scheme suggested 
that six clusters did not significantly deviate from the tree we obtained (1,455 log likelihood 
units, p < 0.55). 

Table 3 shows the Kimura 2P genetic distances among all of the Neotropical Tapirus 
taxa. For all three genes, “T. kabomani” yielded lower genetic distances with regard to 7. 
terrestris than did 7. pinchaque (Cyt-b: 0.018 + 0.003 vs. 0.024 + 0.003; CGI: 0.005 + 0.002 
vs. 0.010 + 0.004; COII: 0.013 + 0.003 vs. 0.017 + 0.004, respectively). The genetic distances 
between “7. kabomani” and 7. terrestris, typical of different populations within a species, are 
not surprising results. However, the small genetic distances between what are traditionally 
considered full species (T.terrestris and T. pinchaque ) are (proof 5). 

Table 3. Kimura 2P genetic distance (Kimura, 1980), in percentages, with standard 
deviations among Neotropical Tapirus taxa pairs for three mitochondrial genes ( COl, 

COII and Cyt-b) 


COI 





Taxa 

T. terrestris 

“T. kabomani” 

T. pinchaque 

T. bairdii 

T. terrestris 

- 




“T. kabomani” 

0.5 % + 0.2 

- 



T. pinchaque 

1.0 % + 0.4 

1.2 % ± 0.4 

- 


T. bairdii 

5.8 %± 1.0 

6.0%+1.0 

5.9% + 1.0 

- 

COII 





Taxa 

T. terrestris 

“T. kabomani” 

T. pinchaque 

T. bairdii 

T. terrestris 

- 




“T. kabomani” 

1.3 % + 0.3 

- 



T. pinchaque 

1.7 % + 0.4 

1.2 % + 0.4 

- 


T. bairdii 

7.8 %+ 1.1 

7.7 %+ 1.0 

8.0 % + 1.1 

- 

Cyt-b 





Taxa 

T. terrestris 

“T. kabomani” 

T. pinchaque 

T. bairdii 

T. terrestris 

- 




“T. kabomani” 

1.8 % ± 0.3 

- 



T. pinchaque 

2.4 % + 0.3 

2.3 % + 0.4 

- 


T. bairdii 

12.5%+ 1.2 

12.1 % + 1.2 

12.4 %+ 1.1 

- 


Figure 4 shows the ML tree for the three concatenated mitochondrial genes with the “71 
kabomani ” sequences (both Brazilian individuals reported by Cozzuol et al., 2013 and three 
specimens sampled by us). Clearly, our tree with 93 specimens showed reciprocal monophyly 
between T. terrestris and 71 pinchaque. Moreover, the alleged “71 kabomani ” is, 
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mitochondrially speaking, a clade within T. terrestris. Thus, our results do not provide 
positive data in favor of T. kabomani as a new and full tapir species (proof 6). The ML tree 
with only the barcoding gene COI employed for species discrimination, also showed “T. 
kabomani ” as a clade within T. terrestris (Figure 5; proof 7). This demonstrates reciprocal 
monophyly between T. terrestris and T. pinchaque. 
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Figure 4. Maximum likelihood (ML) tree for 93 Neotropical tapir specimens (including T. pinchaque, 
T. terrestris, T. bairdii and the alleged new species, “71 kabomani”, reported by Cozzuol et al., 2013) 
for three concatenated mitochondrial genes ( Cyt-b + COI + COIF). T. kabomani was within the T. 
terrestris clade and therefore these mitochondrial genes did not provide positive evidence for T. 
kabomani as a full species. 
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Figure 5. Maximum likelihood (ML) tree for different individuals of the three Neotropical tapir species 
recognized and the alleged new species reported as “71 kabomani” by Cozzuol et al., (2013) at the COI 
gene. “71 kabomani ” was a clade within T. terrestris. 
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The MNJ procedure including the four possible Latin American tapir taxa for the three 
concatenated genes (Figure 6) showed that the haplotype sets of T. bairdii, T. terrestris and T. 
pinchaque are well defined. However, the “T. kabomani” haplotypes seem to be an extreme 
extension from those of T. terrestris (proof 8). 



Figure 6. Median Joining Network (MJN) with the mitochondrial haplotypes found for Tapirus bairdii 
(green), T. terrestris (yellow), “T. kabomani” (brown) and T. pinchaque (blue) with the concatenated 
sequences of three genes (Cyt-b + COI + COII). Red circles are intermediate haplotypes not found. “T. 
kabomani ” was a prolongation of T. terrestris. 

An individual of T. terrestris from the Chiquitania in Bolivia (this was not reflected in 
any phylogenetic tree) connected the “T. kabomani ” and T. pinchaque haplotypes. Some 
temporal split estimations were estimated by this procedure and they are interesting and cited 
here. For example, the temporal split between T. bairdii and T. terrestris and T. bairdii and 
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“T. kabomani ” were 7,960,000 + 69,282 YA (p = 39.8 ± 0.346) and 9,560,000 ± 144,222 YA 
(p = 47.8 + 0.721), respectively. Therefore, on average, the time split between T. bairdii and 
T. terrestris + “I kabomani” was around 8.76 MYA. The temporal split between T. terrestris 
and T. pinchaque was around 3,420,000 + 509,117 YA (p = 17.1 + 2.545) and the haplotype 
diversifications within both species were 3,300,000 + 141,421 YA (p = 16.5 + 0.707) and 
1,208,696 + 194,635 YA (p = 6.04 + 0.973), respectively. In consideration of these data, the 
mitochondrial diversification of T. terrestris preceded T. pinchaque and/or being the 
subsistence of the original haplotypes in T. terrestris higher than in T. pinchaque. This last 
possibility agrees quite well with the fact that the gene diversity is lower in T. pinchaque than 
in T. terrestris. This is probably more due to the greater action in gene drift associated with 
considerably small population sizes in the first species rather than in the second species. The 
temporal split between T. terrestris and “T. kabomani ” was around 1,300,000 + 360,555 YA 
(p = 6.5 + 1.802), considerably lower than that between T. terrestris and T. pinchaque , 
ratifying “T. kabomanr as a clade within T. terrestris (proof 9). 

The MJN applied only to T. terrestris at the Cyt-b gene showed the following picture 
(Figure 7). There were six well-defined haplogroups. Amazon 0 (which corresponds to “T. 
kabomani ” in the other analyses) and Amazon I haplogroups were the most differentiated 
ones. The average temporal split among the six haplogroups was around 1.6 + 0.34 MYA. 
The beginnings of the temporal diversification within each haplogroup were as follows: 0.59 
MYA for the Amazon I, 0.42 MYA for Amazon II and Amazon III, 0.55 MYA for the North 
and 0.33 MYA for the South. 
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Figure 7. Median Joining Network (MJN) with the 80 haplotypes found at the Cyt-b gene for the 141 
lowland tapirs analyzed. Six different mitochondrial haplogroups were found (Amazon 0,1, II. Ill, 
North and South). The Amazon 0 haplogroup corresponds to “71 kabomani”. Red circles indicate 
missing intermediate haplotypes. 


Demographic Changes in Latin American Tapirs 

The mismatch distribution for T. pinchaque (Figure 8) showed no clear demographic 
changes of the total sample or in the Colombian population. But, it did show a significant 
curve related to a population expansion for the Ecuadorian sample ( rg = 0.006, p = 0.011). 
For the Ecuadorian population, this expansion began around 2,900 years ago and the initial 
population size could have been around 4,770-6,691 females. Four of the five demographic 
change tests (Tajima’D, Fu & Li D* and F*, Fu’s F) were significant for the total sample. 
However, the Ramos-Onsins & Rozas R2 test was not significant for any of the three samples 
(Table 4). 

Table 4. Demographic statistics applied to the overall Tapirus pinchaque sample studied, 
to the Colombian sample and to the Ecuadorian sample. + P < 0.05; * P < 0.01, 
significant population expansions 



Tajima D 

Fu & Li 

D* 

Fu & Li F* 

Fu’s Fs 

raggedness 

rg 

R2 

Overall Tapirus 

pinchaque 

sample 

P[D < - 
2.502] = 
0.0000* 

P[D* < - 
4.49] = 
0.0008* 

P[F* < - 
4.49] = 
0.0006* 

P[Fs < - 
7.02 ] = 
0.0032* 

P[rg < 
0.0251] = 
0.1213 

P[R2 < 
0.0814] = 
0.2091 

Colombian T. 

pinchaque 

sample 

P[D < - 
2.279] = 
0.0013* 

P[D* < - 
2.337] = 
0.0263* 

P[F* < - 
2.551] = 
0.0240* 

P[Fs < - 
6.126] = 
0.0074* 

P[rg< 
0.0363] = 
0.2276 

P[R2 < 
0.0923] = 
0.0917 

Ecuadorian T. 

pinchaque 

sample 

P[D < - 
1.825] = 
0.0144* 

P[D* < - 
2.337] = 
0.0253* 

P[F* < - 
2.551] = 
0.0227* 

P[Fs < - 
6.126] = 
0.0029* 

P[rg< 
0.0061] = 
0.0111* 

P[R2< 
0.0944] = 
0.1517 


Total overall T. pinchaque sample (constant population) 
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Figure 8. (Continued). 


Colombian T. pinchaque sample (constant population) 
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Figure 8. Historical demographic analyses by means of the mismatch distribution procedure (pairwise 
sequence differences) for the mitochondrial DNA studied in Tcipirus pinchaque. These analyses were 
applied to: Overall Tapirus pinchaque sample (A); Colombian mountain tapir sample (B); Ecuadorian 
mountain tapir sample (C). 




Table 5. Demographic statistics applied to the overall Tapirus terrestris sample studied, to the haplogroups. 
to the population separated by the Amazon River and to “T. kabomanp . * P < 0.05; ** P < 0.01, 

significant population expansions 



Tajima D 

Fu & Li D* 

Fu & Li F* 

Fu’s Fs 

raggedness rg 

R2 

Overall Tapirus terrestris 
sample 

P[D < -1.642] 

= 0.018* 

P[D* < -4.00] 

= 0.004** 

P[F* < -3.59] 

= 0.004** 

P[Fs < -34.02 ] 

= 0.000** 

P[rg < 0.0035] = 
0.0009** 

P[R2 < 0.0453] = 
0.031* 

Haplogroups 







Amazon I 

P[D < -0.909] 

= 0.176 

P[D* < - 
1.563] =0.089 

P[F* < -1.611] 

= 0.101 

P[Fs<-4.517] 

= 0.008** 

P[rg < 0.051] = 
0.145 

P[R2 < 0.089] = 
0.003** 








Amazon II 

P[D < -1.639] 

= 0.046* 

P[D* < - 
3.307] = 

0.002** 

P[F* < -3.244] 

= 0.003** 

P[Fs < -9.990] 

= 0.000** 

P[rg < 0.051] = 
0.346 

P[R2 < 0.088] = 
0.256 

Amazon III 

P[D < -1.818] 

= 0.015* 

P[D* < - 
3.152] = 

0.005** 

P[F* < -3.198] 

= 0.013* 

P[Fs < -9.691] 

= 0.001** 

P[rg < 0.050] = 
0.360 

P[R2 < 0.089] = 
0.235 

North 

P[D < -1.738] 

= 0.028* 

P[D* < - 
3.242] = 

0.006** 

P[F* < -3.242] 

= 0.010* 

P[Fs < -6.728] 

= 0.008** 

P[rg < 0.051] = 
0.422 

P[R2 < 0.091] = 
0.202 

South 

P[D < -1.714] 

= 0.033* 

P[D* < - 
1.687] = 

0.050* 

P[F* < -1.954] 

= 0.041* 

P[Fs < -2.965] 

= 0.047* 

P[rg < 0.060] = 
0.204 

P[R2 < 0.049] = 
0.019* 

Amazon River as a barrier 







North to the Amazon River 

P[D < -1.215] 

= 0.095 

P[D* < - 
2.489] = 

0.020* 

P[F* < -2.345] 

= 0.026* 

P[Fs < - 
30.424] = 

0.000** 

P[rg < 0.006] = 
0.033* 

P[R2 < 0.060] = 
0.112 

Southern to the Amazon River 

P[D < -1.322] 

= 0.042* 

P[D* < - 
2.872] = 

0.013* 

P[F* < -2.769] 

= 0.012* 

P[Fs < -6.695] 

= 0.024* 

P[rg < 0.006] = 
0.003** 

P[R2 < 0.060] = 
0.025* 

“T. kabomani ” 

P[D < 

0.30550] = 
0.6511 

P[D* < 0.311] 

= 0.6404 

P[F* < 

0.3259] = 

0.6397 

P[Fs < 

1.2411] = 

0.6601 

P[rg < 0.01532 
] = 0.0000** 

P[R2 < 0.079401] 

= 0.0024** 
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The mismatch distribution and the six demographic tests carried out showed a clear 
population expansion for the overall T. terrestris sample studied (Figure 9 and Table 5). All 
six tests were highly significant. There was clear evidence of population expansion in 
different haplogroups. For Amazon I, the mismatch distribution and two out of six tests were 
significant. 

This was the haplogroup where there was less evidence of demographic change. For 
Amazon II, III and North, the mismatch distribution and four out of six tests were significant. 
For South, the mismatch distribution and five out of six tests were significant. Thus, we 
determined noteworthy evidence of population expansion within each haplotype lineage as 
well as for the total T. terrestris distribution. When we only considered two geographical 
groups (North and South of the Amazon River), there was also strong evidence of population 
expansion. For the Southern Amazon group, the mismatch distribution and six out of six tests 
were significant, whereas for the Northern Amazon group, the mismatch distribution and four 
out of six tests were significant. 

Although, the sample of “T. kabomanF is modest, we also analyzed possible 
demographic changes in this taxon (Figure 9 and Table 5). The mismatch distribution and two 
out of six demographic tests showed evidence of population expansion such as was observed 
in T. terrestris, especially in the Amazon I haplogroup. 

Thus, “T. kabomani” should be a dynamic demographic extension of T. terrestris, 
characterized by historical population expansions (proof 10). 

Overall T. terrestris sample (population expansion) 



a 


Figure 9. (Continued). 
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Amazon I (population expansion) 



Pairwise Differences 



Figure 9. (Continued). 
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Southern Amazon River bank (population expansion) 


0.08 


0.04 


0 . 

0 10 20 30 

Pairwise Differences 

C 

“T. kabomani” (population expansion) 

0.2 


0.15 


0.1 


0.05 


0 

0 10 20 30 

Pairwise Differences 

d 

Figure 9. Historical demographic analyses by means of the mismatch distribution procedure (pairwise 
sequence differences) for the mitochondrial DNA studied in Tapirus terrestris and “T. kabomani”. 

These analyses were applied to: Total T. terrestris sample (A); Amazon I, Amazon II, Amazon III, 
North, and South (B); Northern bank of the Amazon River and Southern bank of the Amazon River 
(C); “T. kabomani” (D). 

Spatial Structure 

For T. pinchaque, the overall Mantel test between the log (genetic distance) and the log 
(geographic distance) showed r = -0.127 (p < 0.691 from 10,000 randomizations). Therefore, 
there is no measureable spatial genetic structure in T. pinchaque, at least, in the areas sampled 
in Colombia and Ecuador for mitochondrial DNA. Similarly, for T. terrestris, the Mantel test 
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did not detect any significantly positive correlation between the genetic distances and the 
geographical distances among the exemplars sampled at the global level, within different 
haplogroups and within certain geographical regions at the Cyt-b gene. Some correlations 
were significant but negative (for the global tapir population studied, r = -0.169, p = 0.001; 
Amazon I, r = -0.439, p = 0.001; Amazon H, r = -0.213, p = 0.001; Amazon III, r = -0.107, p 
= 0.023; North, r = 0.007, p = 0.439; South, r = -0.525, p = 0.001; total Colombia, r = -0.070, 
p = 0.007; Northern Colombia-Venezuela, r =-0.134, p = 0.001). Therefore, no isolation by 
distance patterns were found in any of the T. terrestris groups analyzed. Nevertheless, for T. 
terrestris, an AIDA analysis for the entire Amazon basin showed the 1 DC (0-350 km) to be 
significantly positive. The 4 DC (1,150-1,350 km) had a significantly negative autocorrelation 
and again a significantly positive autocorrelation at 6 DC (1,500-1,700 km). This correlogram 
showed the typical resemblance of a double circular dine (Sokal and Oden, 1978), which 
implied a complex pattern of colonization of this area by the lowland tapir. 

In contrast, the T. bairdiVs haplotypes showed a very significant spatial genetic trend (r = 
0.819, approximate Mantel t-test = 3.611, p = 0.0002; out of 5000 random permutations, 5000 
were < Z, 0 = Z and 0 > Z; one-tail probability p = 0.0004), an effect of isolation by distance. 


Discussion 

Genetic Conservation of the Latin American Tapirs 

Clearly, T. terrestris had the highest gene diversity level of all tapir species. This finding 
agrees quite well with the fact that the Latin American Tapirus species had the widest 
geographical distribution and therefore the highest potentially effective numbers. Curiously, 
T. pinchaque, a species with a very restrictive geographic distribution and a small population 
size, presented a gene diversity that was more than three times higher than that of T. bairdii. 
The surprise is due to the fact that T. pinchaque, may be on the brink of extinction (Ashley et 
al., 1996) whereas T. bairdii, has historically occupied a distribution from Southern Mexico 
to the Pacific area of Colombia and maybe Ecuador (Tirira, 2011). Thus, although T. 
pinchaque has a very small census population size and a very restrictive geographical 
distribution, within disturbed and fragmented ecosystems, the species is not impoverished 
from a genetics point of view. This is good news from conservation perspective. 

The geographical distribution of T. pinchaque is still not completely known and there is 
disagreement over its distribution within the scientific literature. In Colombia, populations of 
the mountain tapirs are well described in the Central and Eastern Cordilleras (Lizcano et al., 
2002; Cavelier et al., 2010). Lizcano et al., (2002) estimated around 2,500 mountain tapirs in 
Colombia with around 1,874 animals in seven Colombian National Parks. However, Cavelier 
et al., (2010) estimated a minimum of 543 specimens and a maximum of 579 animals in the 
eight Colombian National Parks where they considered that this species still lives. However, 
the presence of the mountain tapir in the Western Colombian Cordillera is disputed. Acosta et 
al., (1996) and Downer (1997) recorded locations in this Cordillera where the species could 
have lived in the past. Sarria (1993) recorded comments of staff working at Farallones de Cali 
National Park who observed T. pinchaque near Alto de Pance and Alto del Hambre (3,730 
masl). Another possible record of the presence of T. pinchaque within the Western Cordillera 
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was made by De Wilde (1994), who wrote about footprints in the Caramanta Cerro between 
the Departments of Antioquia and Risaralda. However, Lizcano et al., (2002) suggested that 
this species does not inhabit the Western Cordillera. Furthermore, Lizcano and Cavelier 
(2000) considered that such records within the Western Cordillera may indicate the Baird's 
tapir ( T. bairdii) and not T. pinchaque. The Baird’s tapir could also live within the montane 
forest of Central America, up to 3,000 masl. However, Arias-Alzate et al., (2010) recently 
reported the existence of a skull of a T. pinchaque in a Medellin museum sampled in 1911 
from Paramo de Frontino in the Antioquia Department of Western Cordillera. This provides 
evidence of T. pinchaque living within the northern area of the Western Cordillera. 

Lizcano et al., (2002) recorded the presence of the mountain tapir within the Central 
Cordillera from Los Nevados National Park (Caldas and Risaralda Departments) to the South. 
They also noted the mountain tapir within the Eastern Cordillera from Paramo de Sumapaz to 
the South. However, some populations of mountain tapir have been recently detected north of 
Paramo de Sumapaz within the Eastern Cordillera, including at Chingaza National Park and 
12 other northern localities. Paramo de Pisba may be one of these sites (Montenegro, 2002). 
The current absence of mountain tapirs from Northern Central and Eastern Cordilleras could 
be the result of hunting activity (Lizcano et al., 2002). The extreme decline of the female 
population we detected for the Colombian T. pinchaque population at the genetic level with a 
Bayesian sky plot analysis (Ruiz-Garcia et al., 2015b) correlated quite well with the fact of an 
extreme reduction in the potential habitat space of the Andean tapirs. Cavelier et al., (2010) 
showed that the current distribution for Andean tapirs covers 31,400 km , compared to 
205,000 km" in the past. Similarly, the trend in distribution size for tapirs in Colombia is 
much lower today (14,385 km") than in the past (74,556 km"). That is, the current value 
corresponds to 19% of the past distribution. The area where the species occurs is actually a 
collection of 35 forest patches. Most of the fragments are small and only four are larger than 
1,000 km" (Lizcano et al., 2002). These authors showed that only five to six fragments have 
the minimum necessary size (826 km") to maintain at least 150 individuals, the estimated 
number to maintain a viable population in the short term. Thus, the situation of T. pinchaque 
is critical in Colombia. 

In Ecuador, the areas with the highest T. pinchaque population sizes are the Cayambe- 
Coca National Park (Northern Eastern Cordillera) and the Sangay National Park (Central 
Eastern Cordillera). There are other areas of Ecuador, such as the Reserva Ecologica El Angel 
(Carchi Province), Cajas National Park-Bosque Protector Mazan (Azuay Province) and the 
Condor Cordillera, where the mountain tapirs have lived but where there are no records of 
them in the last few decades (Downer, 1997; Tirira and Castellanos, 2001). 

Cavelier et al., (2010) estimated a population size of 543-579 individuals for Colombia 
and 983-1,047 individuals for Ecuador. However, the gene diversity of T. pinchaque in 
Colombia is higher than in Ecuador. Dobzhansky (1971) showed that ancestral populations 
retain more elevated gene diversities than do derived populations. This may support that the 
original T. pinchaque population inhabited the Colombian Andean Mountains and later it 
expanded toward the southern mountains (Ecuador today). Our analyses detected the 
population expansion of the Ecuadorian population. However, more recently, the Colombian 
population decreased due to hunting and habitat fragmentation (high population decrease in 
the last 5,000 YA; Ruiz-Garcia et al., 2015b). The Ecuadorian population has been less 
affected by human activity. Thus, the Colombian population could be in a more dangerous 
situation than the Ecuadorian population. Another T. pinchaque population in critical 
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situation should be the Peruvian population. Cavelier et al., (2010) only estimated between 41 
to 43 mountain tapirs in the Peruvian protected Tabaconas Namballe National Sanctuary. 
Downer’s (1997) work suggests the possibility of mountain tapirs inhabiting Eastern Tapal in 
the Cajamarca Province, south of the Jaen Province, and at the Piura Province. The Peruvian 
population needs a molecular population analysis. 

The gene diversity levels of T. terrestris are the highest estimated for the Latin America 
tapirs and this species seems to not be significantly endangered from a genetics point of view. 
Also, the major part of the haplogroups showed historical population expansions. However, 
T. terrestris is a species threatened due to the subsistence hunting carried out by the Indian 
and “colonos” communities in South America. Another threat comes from the deforestation 
of the Neotropical forests due to extensive agriculture and ranching (Constantino et al., 2006). 
In many areas of South America, this species is more targeted than other animals by hunters. 
Here we provide three examples of hunting: 1- Our first example is the case of 22 Indian 
Izoceno-Guarani communities in the “chaquena” region of Bolivia, who intensively hunt 
tapirs (Barrientos and Cuellar, 2003); 2-A second example is the case reported by Romero et 
al., (2013) in the community of Caura (Guarataro), along the Southern Orinoco River in 
Venezuela. There humans consume 40 % of their proteins from hunted animals. Of the 196.6 
g of proteins consumed by a human per week, 77 g originated from wild hunts and of these, 
14 were from tapirs. Out of 184 tons of wild meat consumed in a year, 44 tons were from 
tapirs. This means that in one year this community hunted 274 tapirs. 3- Third, it is clear that 
the tapir provided the most meat in the Peruvian rainforest compared to all other mammals, 
but unfortunately the hunting model used is non-sustainable. Bodmer et al., (1997a) 
calculated the percentage of production taken by hunters in the Pacaya-Samiria National Park 
within the Peruvian Amazon. For example, the hunted red deer ( Mazama amaricana), 
collared peccary ( Pecary tajacu) and white lipped peccary (Tayassu pecari) yielded values of 
6.8 %, 4.3 % and 19.8 %, respectively. The authors considered that a value higher than 20 % 
could be dangerous to the continued subsistence of a population. The value for the lowland 
tapir in that Peruvian area was 166 %, which showed that the probability of subsistence of 
that population was completely un-sustainable. Also, Bodmer et al., (1997b) showed that in 
one year for three areas of the Pacaya-Samiria National Park, 53 Tayassu pecari, 12 Pecari 
tajacu, 13 Mazama Americana and 20 Tapirus terrestris were killed. But, tapir species have 
the highest biomass (4,000 kg) of all mammals, making up 48 % of the ungulate biomass 
extracted and 35 % of the total mammal biomass extracted. Thus, tapirs are critically affected 
by hunting. Furthermore, this species has a long gestation period of around 13 months, only 
producing one offspring by litter (Padilla and Dowler, 1994). Tapirs reach sexual maturity at 
two years of age (Yamini and Schillhorn, 1988), and females only have one offspring every 
three years. Tapirs also have a very low population density (Bodmer and Brooks, 1997; 
Bodmer et al., 2000). Bodmer and Brooks (1997) reported an overall density of 0.4 ind/km" 
for the Peruvian Amazon forests, while Ayala (2003) found a value of 0.5 ind/km" for the 
Bolivian Chaco. Scientists have offered a range of estimations from 0.028 ind/km" in the 
Northern Colombian Amazon (Chamorro-Rengifo and Cubillos-Rodriguez, 2007) to 1.2-1.6 
ind/km" in the Chiquitanian forest of Bolivia (Arispe et al., 2007) and Brazil (Schaller, 1983). 
For this reason, the conservation of tapirs is crucial within Neotropical environments. 

The reduction of the lowland tapir populations is observable because of its markedly 
reduced distribution. One example of this is the putative subspecies T. t. colombianus, that has 
been extirpated from the lowlands of the Caribe region as well as from the Andean valleys of 
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Colombia. This putative subspecies has been isolated in some localities of the Sierra Nevada 
de Santa Marta (Constantino et al., 2006). T. terrestris is classified in the Appendix II of 
CITES and the IUCN (2003) classified it as Vulnerable. Also, the subspecies, T. t. 
colombianus, is considered to be Critically Threatened (IUCN, 2004). In northern Argentina, 
the tapirr distribution range has been reduced by 46 % in the last century by hunting and 
habitat destruction (Chalukian et al., 2013). 

Thoisy et al., (2010) affirmed that the Amazon River acted as a barrier to gene flow in T. 
terrestris. Our results showed that the Amazon River was only a partial barrier for haplotype 
dispersion for T. terrestris. An interesting additional comment related with this is that the 
jaguar corridor initiative proposed by Rabinowitz and Zeller (2010) and Zeller et al., (2013) 
could be very useful for the lowland tapirs too. The jaguar conservation units could be very 
similar and overlap the region with the largest tapir populations (they are the largest carnivore 
and herbivore mammals in South America). The 182 jaguar corridors proposed by these 
authors (2.6 million km") could also be useful in connecting tapir populations especially in 
the northern and southern geographical areas of South America where the habitat destruction 
carried out by humans is more intense. 

The case of T. bairdii seems to be even more dramatic. Although its geographic 
distribution is wide, it has been dramatically reduced and fragmented in the last two centuries 
and today no more than 6,000 individuals are left in the wild (Ashley et al., 1996). Its 
mitochondrial genetic diversity levels were extremely low compared with other Neotropical 
tapirs. This could mean that this species suffered from a bottleneck and/or the gene drift has 
been more intense on this species by natural or human constrictions. Indeed, this species has 
intensely declined in the last century by habitat destruction and hunting and has been extinct 
in El Salvador and in a major fraction of its original distribution range in Colombian and 
probably completely extinct in Ecuador (Bodmer and Brooks, 1997). 


Systematics of the Neotropical Tapirs and the Case of the Alleged “T. 
Kabomani ” 

All the molecular studies on tapirs, with the exception of that of Couzzol et al., (2013) 
and one of the analyses of Thoisy et al., (2010), determined the split between T. terrestris and 
T. pinchaque to have occurred during the last Pliocene or at the beginning of the Pleistocene. 
Ashley et al., (1996), analyzed the COII gene and determined the temporal split between the 
ancestors of T. terrestris and T. pinchaque to be around 3 MYA. The authors concluded that 
this coincided with the disappearance of the Bolivar Trough and the apparition of the modern 
Panamanian isthmus. Norman and Ashley (2000) added more Perissodactyla species and 
other mitochondrial gene ( 12S rRNA) and estimated new temporal splits. The COII and 
12SrRNA genes supported divergence times between T. terrestris and T. pinchaque of 2.5-2.7 
MYA and 1.5-1.6 MYA, respectively. More recently, Ruiz-Garcia et al., (2012), analyzed the 
Cyt-b gene and showed a Bayesian tree where the ancestors of T. terrestris and T. pinchaque 
diverged around 3.8 MYA (95 % High Posterior Density, HPD: 2.1-4.7 MYA). When a p 
statistic was used on a MNJ network, the two most frequent T. terrestris’s haplotypes 
diverged from the main T. pinchaque’ s haplotype around 1.58 + 0.30 MYA and 1.53 + 0.34 
MYA, respectively. Ruiz-Garcia et al., (2015a) showed a Bayesian tree that supported a 
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temporal split of 3.33 MYA between both T. pinchaque and T. terrestris. Also, Ruiz-Garcfa et 
al., (2015b) performed a mitogenomic study of T. pinchaque and determined the temporal 
split between the ancestors of T. terrestris and T. pinchaque to oscillate from 7.42 to 3.27 
MYA (95 % HPD) in a Bayesian tree. In that study also, with the p statistic from the MNJ 
analysis, and with no priors from other paleontological or molecular studies, assuming that T. 
pinchaque was an ancestral taxon with regard to T. terrestris (the Hershkovitz, 1954’s and 
Haffer 1970’s hypotheses), a split occurred between the taxa around 2.47 + 0.42 MYA. In 
contrast, if T. terrestris is older than T. pinchaque, (Ruiz-Garcia et al., 2012’s hypothesis), the 
temporal divergence between them would have occurred around 3.59 + 0.79 MYA. Our 
current work suggests that the split between species occurred 3.4 MYA. 

These temporal estimates absolutely disagree with what Cozzuol et al., (2013) claimed. 
They reported a temporal split within the clade “T. kabomani”-T. pinchaque- T. terrestris 
ranging from 0.65 to 0.29 MYA (95 % HPD, with the split between T. pinchaque and T. 
terrestris around 0.1-0.3 MYA following these authors). In another approach, Thoisy et al., 
(2010), estimated the divergence between T. terrestris and T. pinchaque to have occurred 
around 0.33 MYA. 

Nevertheless, we are more confident that the first estimates are more accurate than the 
second for three main reasons. (1) In the current work we used more individuals of all the 
Latin American tapir taxa and more genes (better gene diversity estimations in each taxa) than 
the Thoisy et al., (2010)’ work. However, we used a mutation rate, in our MJN analysis that 
was close to the second mutation rate used by Thoisy et al., (2010). They obtained a very 
recent split between T. pinchaque and T. terrestris (0.33 MYA). The temporal split estimates 
provided by Ruiz-Garcfa et al., (2015a,b), ranged from 7.4 to 2.5 MYA. They were obtained 
without any prior or constriction made by the researchers on divergence times or mutation 
rates. Our current estimate of 3.4 MYA is within of the quoted range. (2) The extremely 
reduced divergence time put forth by Cozzuol et al., (2013) for the split among “T. 
kabomani”-T. pinchaque- T. terrestris, and which influences the appearance of “T. 
kabomani”, is easily explained by a temporal constriction imposed by these authors in their 
data. They applied a constriction for the mitochondrial haplotype diversification in T. 
terrestris of 0.13 + 0.1MYA because they affirmed that the oldest fossil records of T. 
terrestris date back to the beginning of the fourth Pleistocene glaciation (0.13 MYA). This 
paleontological constriction helps to explain the difference in temporal divergence between 
South American tapir species noted by Cozzuol et al. (2013) relative to other studies. These 
authors claimed that no clear T. terrestris fossil records dating older than 130,000 YA have 
been located, following Tonni (1992). He described a T. terrestris’s mandible from the 
Lujanense age (upper Pleistocene, 130,000 YA) that he collected from the Colon Department 
at the Entre Rios Province (Argentina). Also, Noriega et al., (2004) and Ferrero et al., (2007), 
found fossil fragments of T. cf. terrestris corresponding to the El Palmar Formation at the El 
Boyero locality (upper Pleistocene) at the Entre Rio Province (Argentina). But, this does not 
mean that they don’t exist, simply that they have yet to be found. Recently, Holanda and 
Rincon (2012) reported two tapir remains in Venezuela. One of them was classified as T. 
terrestris from the Zumbador Cave from the upper Pleistocene, while the second fossil from 
El Breal de Orocual was classified as Tapirus sp. It’s unclear if this last fossil was from the 
Pliocene or from the early Pleistocene. For example, if this last remain is classified as T. 
terrestris, the Cozzuol et al., (2013)’s constriction doesn’t make sense and therefore the 
divergence within the clade “T. kabomani”-T. pinchaque- T. terrestris should be considerably 
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older. Indeed, we believe that some of the oldest South American tapir fossil remains 
classified as Tapirus sp could belong to T. terrestris, especially if a population view is 
adopted more than a typological one, which is more typical of paleontologists. Hulbert et al., 
(2009) showed that the discovery of 75 individuals of T. polkensis in the Gray Fossil Site in 
Eastern Tennessee indicated a unique species. However, it had considerable intraspecific 
variation in the development of the sagittal crest, outline shape of the nasals and the number 
and relative strength of lingual cusps on the PI. These authors concluded that if these fossil 
remains had been found in diverse geographical areas, they would have been considered 
different species. Similarly, Perini et al., (2011) demonstrated that the lower molariform teeth 
size and proportions, used by many authors to define different tapir fossil species, are 
unreliable because they have great population variability. Thus, some early and middle 
Pleistocene tapir fossils should be re-assigned and some of them could be un-differentiated 
from T. terrestris. The findings of Perini et al., (2011) don’t support the conclusion that all of 
these Tapirus fossils are T. terrestris, they simply indicate that some are not, but that others 
could be. Indeed, in our population criteria and in a context of anagenesis, or phyletic 
evolution, some fossil of tapirs, like T. rondoniensis, T. cristatellus and, even, T. 
mesopotamicus could be interpreted as T. terrestris individuals with some small 
morphological differentiated traits. These traits would probably be due to the exposure to 
different environmental conditions in each moment and in each Pleistocene refugium as well 
as due to phenotypic plasticity. Another controversial question between genetics and 
paleontological findings in regard to tapirs is the fact that several authors, from a 
paleontological point of view, claimed that the South American tapirs do not form a 
monophyletic group (Holanda and Ferrero, 2013). These authors maintain that T. pinchaque, 
T. terrestris, T. mesopotamicus, T. rondoniensis and T. cristatellus are paraphyletic with 
evolution in situ in South America—deriving from a similar form such as T. webbi from 
North America during the Miocene. They also claimed that another dispersal event occurred 
from South America to North America by means of a form similar to T. cristatellus, which 
gave origin to the derived forms of North America. This contrasts with the opinion of other 
authors such as Hulbert and Wallace (2005) and Hulbert (2010), who maintained that the 
North American forms were more primitive. From a molecular genetics perspective, Ruiz- 
Garcia et al., (2012), and the current work, found that the two extant tapir species in South 
America only belonged to a unique migration wave, whilst T. bairdii belonged to another 
different molecular group, older than the first and probably more related to the fossil tapirs 
from North America. Thus, the molecular results seem to agree better with the second 
paleontological hypothesis than with the first one. 3- Within the Perissodactyla order, there 
are no cases of an appearance of new species in the last 0.1-0.3 MYA. In Equidae, Xu (1996), 
took into account all of the mitochondrial DNA and estimated the initial splitting within 
Equus to have occurred approximately 9 MYA. Tougard et al., (2001), by means of the 12S 
rRNA and Cyt-b genes, determined that time split to be 12.2 + 2.2 MYA. Recently, Orlando et 
al., (2013), sequencing more than 5,000 genes, determined that all contemporary horses, 
zebras and donkeys originated 4.0^4.5 MYA. Furthermore, they determined the divergence 
time between populations of Przewalski’s and domestic horses to be approximately 0.38-0.72 
MYA (different horse subspecies). In Rhinocerotidae, the split between the current species is 
estimated to be around 26 MYA (isoenzymes; Merenlender et al., 1989) and 21.5 MYA (12S 
rRNA and 16S rRNA genes; Morales and Melnick, 1994). The molecular temporal divergence 
between the two Asian rhinoceros genera was around 25.9 + 1.9 MYA (Tougard et al., 2001). 
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Their paleontological estimates oscillated from 16 to 23 MYA (Carroll, 1988). The molecular 
split for the two Asian rhinoceros species of the genus Rhinoceros (R. sondaicus and R. 
unicornis ) was molecularly estimated to have occurred around 11.7 + 1.9 MYA (Tougard et 
al., 2001). Their paleontological split has been estimated to have occurred 1.6-3.3 MYA 
(Carrol, 1988). The two traditional subspecies of the white rhinoceros ( Ceratotherium simum 
simum and C. simum cottoni) have been estimated to diverge around 1 MYA (Groves et al., 
2010). Henceforth, the claim by Cozzuol et al., (2013) of a divergence time of around 0.1-0.3 
MYA between T. terrestris and T. pinchaque does not agree with that determined for other 
current Perissodactyla taxa. 

With regard to the systematics of T. pinchaque, our genetics data make a remarkable 
contribution. Hershkovitz (1954) found a variable trait in the T. pinchaque dentition, as 
manifested by the first upper premolar, which is variable. Ecuadorian specimens show the 
simple condition, with the cinguloid shelf absent, whereas the first Colombian T. pinchaque 
skull examined shows the premolar as another Tapirus species. For this reason, he 
distinguished two subspecies (T. p. pinchaque and T. p. leucogenys ) within the mountain 
tapir. However, our molecular results showed that Colombian and Ecuadorian specimens 
were mixed independently of their geographical origins disagreeing with the fact that they are 
different subspecies. 

With regard to the systematics of T. terrestris, the existence of the six mitochondrial 
haplotype lineages did not agree at all with the four morphological and geographical 
subspecies determined by Cabrera (1961) and Hershkovitz (1954). That is, there was no 
correspondence among the mitochondrial haplogroups and the morphological subspecies. In 
the geographical area of T. t. aenigmaticus, we found the six different haplogroups detected in 
this study. In the wide distribution of T. t. terrestris, four haplogroups were found (Amazon II 
and III, North and South). Thus, there was no correspondence among the subspecies T. t. 
anigmaticus and T. t. terrestris and the mitochondrial lineages herein showed. We tentatively 
negate the validity of these two subspecies. In the territory of T. t. colombianus, two 
haplogroups were detected, Amazonian II and North. Theoretically, the animals sampled in 
the Colombian Departments of Antioquia, Cordoba and Magdalena (Sierra Nevada de Santa 
Marta) and at the Zulia (Maracaibo Lake) in Venezuela belonged to T. t. colombianus. 
However, one animal from the Antioquia Department belonged to the Amazonian II 
haplogroup. Animals from Meta, Vichada and Guainia Departments (Colombian Eastern 
Llanos and Amazon), French Guyana, Madre de Dios River (southern Peruvian Amazon), 
Yavari River (Western Brazilian Amazon) and even animals from Northern Argentina shared 
haplotypes with the North haplogroup. Therefore, there is no linear correspondence between 
T. t. colombianus and the North haplogroup, which questions the reality of this subspecies. 
The case of these animals from Northern Argentina within this haplogroup is extremely 
strange. It could be that animals with the same haplotypes migrated from the upper part of the 
Western Amazon to the northern and southern distribution ranges. And, for this reason, both 
extreme geographical areas (Northern Colombia and Northern Argentina) share mitochondrial 
haplotypes. Another possibility, which must be remotely considered, is the fact that during 
the XIX century and the beginning of the XX century, people imported tapirs from northern 
areas of South America and liberated them to Northern Argentina for game hunting purposes 
(Martinez, personal communication). This could explain the presence of haplotypes of the 
North haplogroup in Northern Argentina. 
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All the animals we sampled in Southern Brazil, Eastern Bolivia and Paraguay, which “a 
priori” belonged to 71 t. spegazzinni by its geographical distribution and morphological 
features, belonged to the South haplogroup. Therefore all the southern animals belonged to a 
unique mitochondrial haplogroup, with the exception of the Argentinian individuals. In this 
case, it could be a linear correspondence between the South haplogroup and 71 t. spegazzinni. 
However, some animals from the Central Brazilian Amazon and from the Colombian 
Amazon also showed haplotypes of this South haplogroup. It is also interesting to note that 
traditionally only 71 t. spegazzinni has been cited in Bolivia (Anderson 1997). In Bolivia, 35 
specimens have been located representing 24 localities (Salazar-Bravo et al. 2003). We found 
many tapirs of Bolivia belonging to the South haplogroup (related with 71 t. spegazzinii), but 
we also determined exemplars of the Mamore River which belonged to the Amazonian II 
lineage typical of more Northern Amazon areas. 

Due to the complex spatial distribution pattern of 71 terrestris garnered through analysis 
of the Cyt-b gene it is impossible to assign tapirs with unknown origin and living in captivity 
to precise geographical locations. However, we can assign them to one of the six haplogroups 
we found. For instance, the five tapirs from the US zoos belonged to three out of the six 
haplogroups (Amazon I, III and South). Also, the exemplar sampled at the Barcelona Zoo 
belonged to the South haplogroup. Moreover, the animal of unknown origin we analyzed was 
assigned to the North haplogroup. It was a unique tooth that had no recorded origin. This 
complex spatial structure could complicate the application of some biological conservation 
concepts, such as ESUs (Evolutionary Significant Units) and MUs (Management Unit) 
(Moritz and Faith, 1998) to 71 terrestris. 

Our mitochondrial data totally disagree with the results of Cozzuol et al., (2013). They 
claimed that “71 kabomani ” was a full species. We provide 10 population genetics proofs 
against this claim. All of our trees, as well as those from Ruiz-Garcia et al., (2015a,b,c) 
showed that “71 kabomani” is a particular lineage within 71 terrestris. In a mitogenomic study 
with 15 genes and more than 250 T. terrestris individuals, Ruiz-Garcia et al. (2015c) showed 
reciprocal monophyly between 71 terrestris and 71 pinchaque. Also, “71 kabomani”, together 
with the Amazon I haplogroup, were the first haplogroups to diverge within T. terrestris. In 
the current work, with three mitonchondrial genes, the genetic distances between T. terrestris 
and “71 kabomani” were always lower than the genetic distances obtained between 71 
terrestris and 71 pinchaque. 

Kartavtsev (2011) analyzed sequences of the COI gene from 20,731 vertebrate and 
invertebrate animal species. He estimated the average distance data for five different groups. 
He obtained 0.89% ± 0.16% for populations within species, 3.78% ± 1.18% for subspecies or 
semispecies, 11.06% ± 0.53% for species within a genus; 16.60% ± 0.69% for species from 
different genera within a family and 20.57% ± 0.40% for species from separate families 
within an order. For this gene, the genetic differentiation between T. terrestris and “71 
kabomani” was only 0.5 %, clearly within the status of populations within species. Ascunce et 
al., (2003) and Collins and Dubach (2000) reported for Primates at the COII gene, an average 
genetic distance around 5.82% among species within a genus and around 2-4 % for 
subspecies. For this gene, the genetic differentiation between T. terrestris and “71 kabomani” 
was only 1.3 %, clearly within the status of populations within species. Bradley and Baker 
(2001) claimed, for the Cyt-b gene, that values < 2% would equal intra-specific variation, 
values between 2% and 11% would merit additional study, and values >11% would be 
indicative of specific recognition. For this marker, the genetic differentiation between T. 



Mitochondrial Gene Diversity of the Mega-Herbivorous Species ... 


35 


terrestris and “T. kabomani ” was 1.8 %, clearly within the status of intra-specific variation. 
The surprising question is the small genetic distances between T. pinchaque and T. terrestris, 
because these taxa have traditionally been considered full species. Thus, considering 
molecular genetic distances, T. pinchaque, should be considered a sub-species of T. terrestris. 
However, we consider T. pinchaque a full species because no natural or captivity 
hybridization has been reported due to reproductive barriers via stasipatric speciation with 
chromosomal changes (White, 1968, 1978). 

The inaccurate result obtained by Cozzuol et al., (2013) was probably due to the fact that 
these authors only analyzed five samples of T. pinchaque at the Cyt-b gene and only one 
sample at the Cyt-b + COI + COII genes. Indeed, they only analyzed three samples at the Cyt- 
b gene, because two T. pinchaque samples were repeated due to that the two samples of T. 
pinchaque that M. Ruiz-Garcia donated to the CozzuoTs team were also donated by another 
scientist ( that shared samples of the same animals with M. Ruiz-Garcia) to B de Thoisy, who 
latter added his T. pinchaque ’ results to those of Cozzuol et al., (2013). The very small T. 
pinchaque sample used by Cozzuol et al., (2013) probably did not represent all of the 
mitochondrial gene diversity of T. pinchaque. This contrasts with the larger T. terrestris 
sample containing animals from a very wide geographical area and thus retaining a major 
fraction of the mitochondrial gene diversity of this last species, and as the genetic distances 
between both T. terrestris and T. pinchaque were very small, by chance the no representative 
mitochondrial gene diversity of T. pinchaque of the small sample of Cozzuol et al., (2013) 
was nested within the gene diversity of T. terrestris. Nevertheless, as soon as we enlarged the 
T. pinchaque sample, collected from a wider sampling area, this phenomena disappeared. 

Cozzuol et al., (2013) also provided alleged evidence of morphological and 
morphometric differences of “T. kabomani ” from the remaining living and fossil South 
American tapir species. In this work, we provide molecular proofs that do not support the 
highly speculative and inconclusive conclusions of Cozzuol et al., (2013) that favor “T. 
kabomani.” We also provide several comments on the alleged morphological differences 
between “T. kabomani” and T. terrestris. In a recent study, Ruiz-Garcla et al., (2015c) 
analyzed several T. terrestris populations throughout Northern Colombia and different areas 
of the Amazon basin in Colombia, Peru, Brazil and Bolivia from a craniometric perspective 
(around 160 skulls). Most of these populations showed significant statistical differences 
regarding the two areas of the skull where Cozzuol et al., (2013) found the greater divergence 
of “T. kabomani” with regard to T. terrestris (the position of the frontal-parietal suture related 
with the beginning of the sagittal crest and frontals broad and more inflated behind the nasals 
in “T. kabomani”). The significant differences of these skull from those studied by Ruiz- 
Garcla et al., (2015c) within diverse T. terrestris populations were related to different 
ontogenetic patterns within each population as well as to the age composition of each 
population. The statistical differences were extreme among some of the seven T. terrestris 
populations studied for different morphometric traits (1- the Napo River area in the Northern 
Peruvian Amazon, 2- Yavari River and neighborhoods in the Amazonian border of Colombia, 
Brazil and Peru, 3- Sierra Nevada de Santa Marta in Northern Colombia, 4- Meta, 5- Caqueta, 
and 6- Vichada departments, all three in Colombia, and 7- the Mamore River in the Bolivian 
Amazon). However, we don’t claim that each one of these populations was a different tapir 
species, such as Cozzuol et al., (2013) claimed with “T. kabomani”. 

Additionally, we show two other morphological proofs which don’t support the claims of 
Cozzuol et al., (2013). For example, supposedly T. kabomani,” with its lower sagittal crest 
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and wider post-nasal exposition of frontals, is much smaller and has a darker color than the 
“typical” T. terrestris. In Figure 1A, we show a tapir captured in the Amacayacu River 
(Colombian Amazon) with a “T. kabomani ” mitochondrial haplotype. We also show, in 
Figure IB, a “typical” morphological T. terrestris individual from Tena, at the upper Napo 
River in Ecuador. But, it had a “71 kabomani ” haplotype. Both animals exhibited different 
morphotypes to that claimed by Cozzuol et al., (2013) for “T. kabomani” . For example, the 
Ecuadorian animal showed a considerable size and its color was brown. In contrast, as we 
show in Figure 10, in our travels to obtain tapir samples, we have found small-sized lowland 
tapirs. Figure 10A shows a pair of small-sized tapirs sampled in the Colombian Amazon, 
whereas Figure 10B shows a small tapir sampled in the Ecuadorian Amazon. We inspected 
the dentition of all three examplars and determined them to be adults because each had at 
least two erupted molar. Nevertheless, no one showed “T. kabomani” haplotypes. The 
Colombian exemplars belonged to the Amazon II and III haplogroups, while the Ecuadorian 
individual belonged to the Amazon I haplogroup. Indeed, these animals were of a size smaller 
than the two individuals in the photo published by Cozzuol et al., (2013) and that were 
defined as “T. kabomani ” Additionally, and mistakenly, they affirmed the existence of photos 
of “T. kabomani” from Colombia. This was a mistake because we provided the authors with 
two Colombian samples that they claimed to be “T. kabomani ” Furthermore, we did not 
provide any photos of these animals. Most likely, they never asked for the morphology of 
these exemplars to verify a link between these “T. kabomani ” haplotypes and the size and 
dark coat color of these individuals. Complementary, we have observed some lowland tapir 
adults with small sizes and certain bones deformations as if they were affected by 
malnutrition or some kind of genetic syndrome (for instance, small head and body but large 
feet). 



Figure 10. (Continued). 
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Figure 10. Some T. terrestris adult exemplars with very small sizes but without “T. kabomani” 
mitochondrial haplotypes. Two small adult individuals from Leticia, Colombian Amazon (A); one 
small adult individual from Macas, Ecuadorian Amazon sampled by the first author (B). 

Dobzhansky (1937), in his seminal book for the Neodarwinism synthesis, showed that in 
evolution the approaches should be of a population nature rather than a typological one. It is 
quite understandable that a paleontologist (such as is the first author of Cozzuol et al., 2013), 
working with the scarce material that provides the fossil record on most of the occasions, 
keeps a typological morphology vision of the species concept. Indeed, many paleontologists 
only accept a cladogenetic mechanism of speciation following the principle of the punctuated 
equilibrium theory (Eldredge and Gould, 1972; Gould and Eldredge, 1993). This means that 
the punctuated equilibrium can become tautological because paleontologists define fossil 
species on the basis of morphological change. So, it is a trivial problem to observe a strong 
correlation between speciation and morphological changes. However, population geneticists 
and evolutionary biologists (including Darwin) know the existence of rapid morphological 
evolution in a hypothetical ancestral form without speciation (phyletic transformation or 
anagenesis). Therefore, the existence of a certain morphological and morphometric skull 
variability in the lowland tapirs did not mean the real existence of different species. 
Furthermore, the “T. kabomani” case is a current living case, and not a paleontological one. 
With a living species, there are many other aspects to be analyzed in addition to changes in 
morphology before concluding on the existence of a new species. And, changes in 
morphology are the only aspects that can be detected in the fossil record. This project and 
other studies show that the samples sizes and the interpretation of mitochondrial and 
craniometric data by Cozzuol et al., (2013) are insufficient to define a new tapir species in the 
Amazon. Additionally and unfortunately, Cozzuol et al., (2013) did not provide sufficient 
geographical (they contrarily argue that “T. kabomani ” is simpatrically extended by all the 
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Amazon with T. terrestris), ecological, ethological or reproductive (post and/or prezygote) 
barrier information to claim to “T. kabomanF as a new full species. This does not mean that 
other tapir species do not exist, but that the provided data are insufficient and misunderstood. 
In what-ever-the-case, nuclear DNA and especially comprehensive karyological studies (due 
to the possibility of stasipatric speciation; White, 1968, 1978) should be carried out before 
claiming the existence of new species of tapirs. 


The Impact of the Pleistocene Changes on the Phylogeography of T. 
pinchaque and T. terrestris 

The temporal split between the ancestors of T. terrestris and T. pinchaque occurred 
around 3.4 MYA. And, the beginning of the T. terrestris mitochondrial diversification was 
around 3.1-3.3 MYA. These times match up well with the end of the Pliocene and the 
beginning of the Pleistocene (2.6-1.8 MYA; Van der Hammen, 1992), which were extremely 
cold. The average temperature descended 4 °C + 2 °C and rainfall declined (500-1000 mm 
less than today) (Van der Hammen, 1992, 2001). The Pleistocene had 25 main glacial- 
interglacial cycles with climatic cycles changing every 20,000 Y for the first million years 
and thereafter every 100,000 Y. These climatic cycles had oscillations in temperature of 10 °C 
at 2,500 m above sea level (masl). This first lowland tapir haplotype split (3.0-3.3 MYA) also 
coincided with the last formational phase of the Central Andes. The entire Andean chain 
between Cajamarca and Huancavelina was formed during this epoch by volcanic activity. The 
lagoons in Huancho (Northern Lima) and the bay in La Ventanilla Beach also formed during 
this time (Feon-Canales, 2007). The diverging process of the Amazon 0 haplogroup or “T. 
kabomanF and the haplotype diversification within T. pinchaque began around 1.3-1.2 MYA. 
This coincides with the Pre-Pastonian glacial period (1.3-0.8 MYA), which was extremely 
dry and cold. While the Buenos Aires province where mammalian fauna ( Categonus, 
Cyonasua, Clyomys, etc) lived was basically tropical from 2-1.3 MYA, this climate later 
changed. This area was Patagonian Steppe-like around 1.3-0.8 MYA which agrees quite well 
with a dry and cold climate (Forasiepi et al., 2007). Examples of mammalian fauna living at 
this time were guanaco, Lestodelpys and Lyncodon. The diversification process within some 
of the T. terrestris haplogroups occurred around 0.5-0.3 MYA which showed that the same 
climatic events affected these lowland tapir haplogroups. This correlates to the beginning of 
the Kansas-Mindel glacial period, 0.4 MYA. Much of the haplotype proliferation we note in 
the haplogroups identified occurred around 0.2-0.06 MYA and in 0.03-0.01 MYA. This can 
be explained by the extreme climatic changes in these epochs. The last glacial-interglacial 
super-cycle began around 0.13-0.15 MYA (Wiirm. Wisconsin, Vistula, Gambliense). From 
130,000 to 80,000 YA (Eemiense period), the climate was hot, the precipitation greater, and 
there were extensive swampy forests of Aliso, Vallea and Weinmannia (Van der Hammen, 
1992). But, an extremely cold climate (upper Pleniglacial period) began around 70,000 YA 
which could have had a fragmentation effect on the lowland tapir haplogroups. From 30,000 
YA to 10,000 YA there were many periods of extreme dry and cold conditions. For example 
from 30,000 YA to 23,000 YA, there was a very cold period in Europe named Dryas I. 
During this period, the Fuquene Fagoon near the current Bogota (Colombia) dried (29,000 
YA; Van der Hammen, 1992) and Mac Neish (1979). This determined changes in the soil 
acidity and the types of pollen at the Pikimachay Cove found within Peru which are related to 
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the extreme cold 23,000 YA. Later, around 19,000-16,000 YA, the Last Glacial Maximum 
(LGM or Dryer II) occurred, which had the most extensive distribution of snow and ice in the 
Central Andes within the last 200,000 YA. Metivier (1998) estimated that about 18,000 YA, 
the glacial extension in North and Central Andes was around 371,306 km", whereas it is 
currently around 3,220 km" (only 1% of the LGM). Also, there were intense cold periods 
from 14,000 YA to 10,000 YA (with alternating hot periods). Rodbell and Seltzer (2000) 
found moraines (the front of a glacial) up to 3,170 masl (12,000 YA) at the Cordillera Blanca 
(San Martin Department, Peru). Comparatively, today the glacial front is 4,600 masl. Also, 
Rodbell (1991) determined that the other areas of the Cordillera Blanca within the limits of 
the Departments of Huanuco and Ancash (17 glaciers) had snow lines at 4,200 masl, whereas 
today this limit is 500-900 m higher. Similarly, Seltzer (1987, 1990) determined the snow 
lines in the Huaytapallama, Junin and Huancavelina Departments (Peru) 10,950 YA to be 
1,400 m lower than today. These climatic conditions could have caused the last lowland tapir 
haplotypes to diversify. 

Several authors (Ab'Saber, 1982; Brown et al., 1974; Fjeldsa, 1994; Haffer, 1969, 1974, 
1982, 1987, 1997, 2008; Prance, 1982, 1996; Prance and Lovejoy, 1985; Terborgh, 1992; 
Van der Hammen, 1975; Vanzolini, 1970, 1973, 1992; Vanzolini and Williams, 1970; 
Whitmore and Prance, 1987) claimed that these climatic events and subsequent dry/humid 
cycles created in the Amazon basin are the result of Milankovitch cycles. During the dry 
periods rainforest communities split into isolated refugia separated by savannah or arid 
pampas. This hypothesis was originally established for the Pleistocene, but later expanded to 
include the Miocene-Pliocene periods as well and it was named the Refugia Hypothesis. This 
hypothesis could explain the appearances of the different haplogroups within T. terrestris. 
During each dry period the haplogroups were formed, but later during the humid periods the 
tapirs migrated allowing its different haplogroups to diversify in different geographic areas. 
Along with this, the relationship between the original refugia areas was lost which helped 
create the geographical distribution of tapirs we see today. However, some Pleistocene 
refugia could have played an important role in the origins of the haplogroups. For example, 
the Amazon I and the “T. kabomani” haplogroups could have radiated from the Napo refugia 
(Northwestern Amazon), while the Amazon II and III haplogroups could have radiated from 
the Napo and Inambari refugia (Southwestern Amazon). The South haplogroup could have 
expanded from the Inambari or from the Rondonia refugia (between the Madeira and the 
Tapajos rivers), whilst the North haplogroup could have originated from the Bolivar refugia 
within the Guiana area. Another hypothesis, called the Recent Lake hypothesis (Klammer, 
1984; Marroig and Cerqueira, 1997; Nores, 1999, 2004; Sombroek, 1966), could also support 
the formation of these lowland tapir haplogroups. This hypothesis claims that most of 
Amazonia was covered by a huge lake or lagoon during the Pliocene, and successively 
smaller portions of Amazonia were covered during a series of assumedly high sea level events 
during the Pleistocene. These created different islands and archipelagos within this Amazon 
lake. The different haplogroups could have been created on these islands. However, like the 
tapirs, who are exceptionally good swimmers, these haplogroups could have migrated 
throughout aquatic systems and expanded their original geographical areas. Therefore, both 
the Refugia and the Recent Lake hypotheses could support the existence of lowland tapir 
haplogroups. The significant and complex spatial structure found for the lowland tapirs in the 
Amazon basin by means of AIDA agrees quite well with both of these hypotheses. 
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Using Bayesian sky plot analyses, we determined the population changes over the last 
few thousand years for some of the lowland tapir haplogroups and geographical areas studied. 
Amazon I suffered a population decline 5,500 YA coinciding with one of the dry periods of 
the Holocene after the Optimum Climaticum (OP). Following Rothlisberger (1987), around 
6,300 YA, there was a significant increase in temperature especially in Southern South 
America as well as in the Central Andes, evidenced by Ois levels in Huascaran snow. This 
dry period around 5,500 YA was detected in the Amazon, Caqueta and lower Magdalena 
River basins as well as in Andean lagoons in Colombia and Peru (Thompson et al., 1995; Van 
der Hammen & Cleef, 1992; Van der Hammen, 2001). However, in the last 1,400 years, this 
haplogroup again increased. Amazon II and III revealed a population decrease 75,000 YA 
coinciding with the ending of the Eemiense inter-glacial period and the beginning of the 
upper Pleniglacial period (Wiirm I; Van der Hammen, 1992). These two Amazon haplogroups 
began to again expand around 12,000 and 14,000 YA respectively. This could have been the 
result of finishing the coldest and driest period of the fourth Pleistocene glacial period, a 
period that is also related to the last major extinction period for mammals. 

In contrast, North decreased 14,000 YA, agreeing with the massive extinction of 
mammals across the Earth including in South America. This corresponds with the Younger 
Dryas (Dryas III), typical of Northern Europe and Scandinavia (Clapperton, 1993). This 
means that the end of this extremely cold and dry period was not the same for all of South 
America and therefore the diverse haplogroups of tapirs were differentially affected. The 
Amazon’s climate was not very affected by these drastic conditions—leading to an increase 
in some Amazon haplogrpups. However, in Northern South America, the Younger Dryas was 
more drastic and the North haplogroup was negatively affected. Nevertheless, this lineage 
increased around 3,000 YA, just when the average temperature reached a value similar to 
what it is today (Van der Hammer, 1992). Finally, South suffered a population declination 
around 6,000 YA coinciding with the drier period between 7,000-5,500 YA that we 
commented on above. 

When the data were analyzed by geographical area, only one region showed an important 
population declination in the last 5,000 years. This was the lowland tapir population of 
Northern Colombia and Northwestern Venezuela. This strong population decrease, detected 
by genetic methods, agrees quite well with the situation of this population, which was 
considered Critically Threatened in 2004 by the IUCN (Constantitno et al., 2006). 

In contrast to what was expected, in the Andes Mountains, the Pleistocene climatic 
changes did not create a phylogeographic structure in T. pinchaque. This is a conundrum 
because in other Andean mammals, such as the Pampas cat ( Leopardus pajeros) (Cosslos et 
al., 2010; Ruiz-Garcia et al., 2013), the Andean cat (L. jacobita) (Cosslos et al., 2012; Ruiz- 
Garcla et al., 2013) or the spectacled bear (Tremarctos ornatus) (Ruiz-Garcia, 2003, 2013; 
Ruiz-Garcia et al., 2003, 2005), the phylogeographic structure for different kinds of molecular 
markers were significantly marked. This may support that the haplotype diversification within 
the current T. pinchaque occurred more recently than it did in T. terrestris. It’s also possible 
that the gene flow capacity of T. pinchaque was relatively higher for this species in the 
Andean cordillera than for T. terrestris in the Amazon lowlands. The population expansion 
detected in the Ecuadorian T. pinchaque population could correlate with this last idea. If the 
mitochondrial diversification of T. pinchaque occurred much more recently than in T. 
terrestris, this could agree with an event of anagenesis, or phyletic evolution, with the 
original ancestor more related to T. terrestris than to T. pinchaque. T. pinchaque could have 
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appeared via adaptation and exploitation of a new habitat within the Andean Highlands from 
a lineage, or population, of T. terrestris. It could have also appeared via stasipatric speciation. 
However, all of our phylogenetic trees do not agree with this perspective. 

To complement these new mitochondrial results future studies should include sequences 
of HLA markers, autosomic and sexual chromosome introns and other nuclear DNA genes. 
Also, future studies focus on sequencing of tapir fossil remains may be helpful in our 
understanding of the evolution of the current mega-herbivores of South America. 
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